Finite-temperature dynamics of matter-wave dark solitons in linear and periodic 
potentials: an example of an anti-damped Josephson junction 
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We study matter-wave dark solitons in atomic Bose-Einstein condensates at finite temperatures, 
under the effect of linear and periodic potentials. Our model, namely a dissipative Gross-Pitaevskii 
equation, is treated analytically by means of dark soliton perturbation theory, which results in a 
Newtonian equation of motion for the dark soliton center. This reduced model, which incorporates 
an effective washboard potential and an anti-damping term, constitutes an example of an anti- 
damped Josephson junction. We present a qualitative (local and global) analysis of the equation of 
motion. For sufficiently small wavenumbers of the periodic potential and weak linear potentials, the 
results are found to be in good agreement with pertinent ones obtained via a Bogoliubov-de Gennes 
analysis and direct numerical simulations. 
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I. INTRODUCTION 



During the last years, matter-wave dark solitons have 
been the focus of many research efforts in the physics 
of ultracold atoms, due to their ubiquitous presence in 
the formation and nonlinear dynamics of atomic Bose- 
Einstein condensates (BECs) 0, • Importantly, these 
efforts have been driven by considerable experimental 
developments along this direction. In particular, ex- 
perimental investigations on matter-wave dark solitons 
started over one decade ago [3- 6]; however, since the soli- 
ton lifetime was rather limited (either due to the instabil- 
ity of dark solitons in higher-dimensional setups, or due 
to the detrimental effect of strong thermal fluctuations), 
these early experiments were not able to characterize 
many of the soliton properties. Nevertheless, there has 
been a new generation of more controllable experiments 
which have been carried out at extremely low 
temperatures and -at least in some cases- in more quasi 
one-dimensional (ID) setups, so as to enable the defini- 
tive observation of robust moving, oscillating and poten- 
tially interacting dark solitons, in very good agreement 
with the corresponding theoretical predictions. Further- 
more, dark solitons were also observed recently in multi- 
component BECs, where they were coupled either with 
dark or bright solitons, thus forming "dark-bright" 
Ho| or "dark-dark" [3, Ell soliton states, respectively. 

In the same context, and since BEC experiments are 
obviously performed in the presence of thermal fluctua- 
tions, the study of matter- wave dark solitons in finite- 
temperature BECs is a quite relevant and interesting 
problem. Pertinent theoretical studies have been per- 
formed in various settings and using different approaches 
(see, e.g., Refs. [l7l - [25j ). A relevant - and analytically 
tractable - mean-field model that has recently gained at- 
tention [l8|, US [H|-[25| ( a l so i n the context of vortices 
(26j) is the so-called dissipative Gross-Pitaevskii equation 
(DGPE). This model, which was first introduced phe- 



nomenologically [27j and later was justified from a mi- 
croscopic perspective [28[ , can describe accurately finite- 
temperature-induced soliton decay: results stemming 
from a pcrturbative study of soliton dynamics in the 
framework of DGPE, compares favorably to ones ob- 
tained by the more accurate stochastic Gross-Pitaevskii 
model Ho, HH, [23| (the latter, is a grand canonical theory 
of thermal BECs, containing damping and noise terms 
which describe interactions of low-energy atoms with a 
high-energy thermal reservoir (28j). 

Our aim in the present work is in a sense two-fold. On 
the one hand, we aim to study matter-wave dark soli- 
tons in thermal BECs in an experimentally relevant and 
realizable setting. On the other hand, we are aiming 
to reverse engineer an interesting dynamical system that 
has been studied in the context of damping (and external 
constant driving) in some detail in the context of Joseph- 
son junctions; see relevant details in Ch. 2 of Ref. [29[ 
and Ch. 8 of Ref. [lo] for physical and mathematical as- 
pects, respectively. The present setting of BECs enables 
the possibility of analyzing this dynamical system in the 
presence of an anti-damping mechanism induced by the 
coupling of the dark soliton to thermal fluctuations. 

In particular, we focus on the examination of matter- 
wave dark soliton dynamics in the presence of the follow- 
ing effects: (i) coupling to thermal excitations, described 
in the framework of the DGPE model discussed above; 
(ii) a (weak) periodic, so-called optical lattice (OL) [3lj |. 
potential [321 ] ; (iii) an (also weak) electric or gravita- 
tional field, inducing a linear external potential - see, 
e.g., Ref. [33f for an experimental realization. 

The steps of our analysis, and the structure of our pre- 
sentation, are as follows. First, starting from a DGPE 
which encompasses the above potentials (linear andperi- 
odic), we employ dark soliton perturbation theory [2|,|3J- 
|36[, to derive an equation of motion for the ensuing dy- 
namics of the soliton center position (section II). Moti- 
vated by its direct analogy with the equation describing 
Josephson junction dynamics (which has exciting fea- 
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tures such multi-stability and hysteresis), we will offer 
a systematic local and global analysis (fixed points and 
phase plane structure) of the pertinent dynamical sys- 
tem. Then (in section III), we will present results of sys- 
tematic numerical simulations, comparing our analytical 
results with ones obtained from direct numerical integra- 
tion of the DGPE. Finally (in section IV), we will sum- 
marize our results and present our conclusions, as well as 
some interesting possibilities for future extensions of this 
work. 



II. THEORETICAL APPROACH 

A. Perturbation theory for dark solitons. 

The dynamics of a quasi-lD condensate (elongated 
along the ir-direction) under the presence of thermal fluc- 
tuations, can be described by the following dimensionless 
DGPE [Hill : 



(i - j)u t + -u xx 



\u\ 2 u + fxu- V(x)u = 0. (1) 



Here, subscripts denote partial derivatives, u(x, t) is the 
condensate wavefunction, (i is the chemical potential, 
and the term —jut accounts for finite-temperature ef- 
fects (note that the parameter 7 scales with temperature 
T according to a power law, 7 cx T a , with 1 < a < 4) 
[20I [22l . [23| . Finally, the external potential V(x) is as- 
sumed to of the form: 



V(x) = -Ex + Vo cos(2fcr), 



(2) 



where, in the right-hand side, the first term represents an 
electric or gravitational field of strength E, and the sec- 
ond term describes an optical lattice potential of strength 
Vo and wavenumber k [33j. 

Without the external potential and the thermal term, 
i.e., when V{x) = 0, 7 = 0, Eq. ((T|) possesses a dark soli- 
ton solution on top of a constant background of density 
/1; this solution is of the form 



u(x,t) = y^[cos</>tarih(Jf) 



(3) 



where X — ^fJL cos <f>[x — Xo(t)] and xo(t) — v / /I(sin(/))t 
is the soliton center. The amplitude and velocity of the 
soliton are respectively given by ^f\x cos <p and */jusin0, 
while <j> (with \<p\ < ir/2) is the soliton phase angle (so- 
called "black" and "grey" solitons correspond to <fi = 
and <j) 7^ 0, respectivel y). 

Following Refs. [20|, |22|, we now seek a solution of Eq. 
([T]) in the form u(x,t) = Ub(x) exp[— i6(t)]v(x, t), where 
Ub(x) and 0(t) denote the background amplitude and 
phase, respectively, while the unknown complex function 
v(x,t) denotes a dark soliton. We assume that the con- 
densate dynamics involves a fast relaxation scale to the 
ground state, which can be approximated in the frame- 
work of the Thomas-Fermi (TF) approximation, as fol- 
lows: Ub w •\/ max {/ i — V( x )t®}- Then, the evolution of 



the dark soliton on top of this ground state, is described 
by the following perturbed nonlinear Schrodinger (NLS) 
equation: 



1 

w t + -v xx 



(\v\ 2 - l)v = P[v], 



(4) 



where we have used the scale transformations t — > /it, 
x — > y/Jix, and V — > V//J,, and the functional perturba- 
tion P[v] is given by: 



1, 



P[v] = (1 - \v\ 2 )vV + -V x v x + jv t 



(5) 



Below we will treat P[v] as a small perturbation (at 
least in the spatial region where the dark soliton dy- 
namics takes place), assuming that the thermal fluctua- 
tions (governed by the parameter 7) and the strengths Vo 
and E of the linear and OL potentials are weak, namely, 
7 1, Vo <?; 1 and £« 1. Under these assumptions, 
we can then apply to Eq. the perturbation theory for 
dark solitons, originally developed in Ref. |34j and then 
adapted to dark solitons in BECs in Ref. |35l ] (see also 
Ref. [36j], for an application to BECs confined in OLs, 
and Ref. 0] for a review). 

The fundamental premise underlying the dark soliton 
perturbation theory is that, to the leading order of ap- 
proximation, the soliton evolves adiabatically; in other 
words, it preserves its form of Eq. ([3]), yet its parameters, 
namely the phase angle <f> and the position of the center 
xq become slowly-varying functions of time. More specif- 
ically, the solution of Eq. (Q} is taken to be of the form 
v = cos ip(t) tanh[cos ip(t)(x — Xo(t))] + ishnp(t), where 
x o{t) = J sin<p(i)di and tp(t) denote the time-dependent 
center and phase angle of the dark soliton, respectively. 
These functions are found to evolve according to the fol- 
lowing equations 0, l34Tj3l| : 



1 



2 cos 2 ip sin (p 
x — sm (p(t), 



Re 



+ OO 



P 



(6) 
(7) 



where overdots denote time derivatives and asterisk de- 
notes complex conjugate. Evaluating the integral in 
Eq. ([B]), and considering the case of almost black soli- 
tons with sufficiently small soliton phase angles, we ne- 
glect terms of order 0(d<p/dt). Then, we can derive from 
Eqs. dU)-© an equation for the soliton center; this equa- 
tion, when expressed in the original variables, has the 
following form: 



2 1 

-7/xio + -E + A k kV a sm(2kx a ), (8) 



where 



Ai 



nk 



Tik 



v /' V 1 + M ) CSCh V V /' 



(9) 



and Af. w 1 for small k. Note that Eq. {5} has the form of 
a Newtonian equation of motion of a classical unit-mass 
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particle, in the presence of (anti)damping, of strength oc 
(2/3)7^, and in the presence of the effective "washboard" 
potential: 



U(x) = i (A k V cos(2/cx) - Ex) 



(10) 



Alternatively, for the purposes of our analytical consid- 
erations, Eq. © can be expressed as the following dy- 
namical system: 



y, 

2 1 



AkkVo sm(2kx). 



(11) 



It is quite remarkable that this dynamical system 
draws a direct analogy between the dynamics of a dark 
soliton in a thermal BEC under the combined effect of 
gravity and an optical lattice, and that of a Josephson 
junction. In the latter, the dynamical variable is the rel- 
ative phase between the two superconducting elements, 
whose time derivative is associated with the voltage in the 
junction. The bias current is analogous to our linear field 
here, while the Josephson current is the one that provides 
the periodic potential. Lastly, the lossy resistor part in 
the original Josephson junction setting (i.e., 7 < 0) is 
replaced here by the anti-damping effect (i.e., 7 > 0) 
induced by thermal fluctuations. It is also quite note- 
worthy that the original experimental setting of Ref. [H[ 
(see also Ref. [37j), which is the closest to our present 
proposal, was precisely designed to examine macroscopic 
quantum interference between the atomic packets in ad- 
jacent wells of the lattice in close analogy to the so-called 
alternating-current (ac) Josephson effect. Here, we take 
this analogy one step further by parallelizing the motion 
of a dark soliton through this potential environment to 
the evolution of the phase of the Josephson junction. 



B. Qualitative analysis of the dynamics 

Let us now analyze in detail the dynamics of the system 
of Eqs. (fTTj) . First, we note that due to the periodicity 
of the nonlinearity in the system dill) , we will restrict 
the analysis to the interval < x < n/k (one period of 
the periodic function within effective potential). Second, 
it is relevant to observe that two different possibilities 
occur: the existence and non-existence of fixed points of 
Eqs. (fTTj) . Obviously, fixed points exist for E/2AkkVo < 
1, i.e., for E/2AkkVo < 1 (existence of two fixed points) 
or E/2AkkVo = 1 (existence of one fixed point). First, we 
will provide the local analysis for the case E/2AkkVo < 1. 
This case, as well as the one with E/2AkkVo — 1, and the 
scenario where no fixed points exist, will be investigated 
in the next section devoted to the global analysis of the 
dynamics. 



1. Local analysis 

In the case E/2AkkVo < 1, there exist two fixed points: 
x = (x-,0) and x + = (x + , 0). The x-coordinates of the 
fixed points x± satisfy the transcendental equation 



sin(2kx±) 



E 



2A k W ' 



(12) 



and correspond to the local minimum and the local maxi- 
mum of the potential energy U{x) respectively. To study 
the stability of the two fixed points, we first consider the 
linearization of the system (|lip around them, whereby 
the Jacobian matrix reads: 



(13) 



where a = k(4A 2 k k 2 V^ - E 2 ) 1 / 2 . The fixed point x + 
[where its cc-coordinate, 5+, is the maximum of the po- 
tential energy U(x)] has eigenvalues: 





( 






^ ±a 





Ai5(7) = 57M±y~7V+a. 



(14) 



Since a > 0, it is clear that > and X^' < and, 
hence, x+ is always unstable: in fact, it is a saddle point, 
independently of whether 7 = (corresponding to the 
Hamiltonian variant of the model) or 7 > (correspond- 
ing to the finite-temperature case). 

On the other hand, the fixed point x_ [where its x- 
coordinate, x + , is the minimum of the potential energy 
U(x)] has eigenvalues: 



k (+) 



1 



,2, ,2 



9 



(15) 



In the Hamiltonian case 7 = (where the Hamiltonian 
energy T-L(x,y) = \y 2 + U(x) is conserved), the above 

eigenvalues are purely imaginary, i.e., \\ <] = Ha, and 
thus x_ is a stable center, in this case, the dark soli- 
ton center performs oscillatory motions around the min- 
imum of the effective potential. However, in the finite- 
temperature case of 7 > (where the Hamiltonian energy 
%{x, y) is not conserved), we need to separate two cases. 
If ^7 2 M 2 < ct (i.e., for sufficiently low temperatures), the 
eigenvalues are complex with positive real part and, thus, 
the fixed point x_ is an unstable spiral. Physically, this 
means that the dark soliton center will perform a motion 
with an oscillatory growing amplitude around the mini- 
mum of the effective potential, with a growth rate ^7/i 
and an oscillation frequency: 



WcfF 



a - ^7 V 



(16) 



On the other hand, if 7 is sufficiently large (such that 
^7 2 /z 2 > a), the fixed point x_ becomes an unstable 
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node. Physically, this means that the dark soliton cen- 
ter will perform a motion with an exponentially grow- 
ing amplitude. Note that, in both cases, the soliton will 
eventually decay at the rims of the condensate. 

At this point, we should mention that the local anal- 
ysis highlights the emergence of, at least, two different 
bifurcations. First, the transition from 7 = to finite 
7 > is a bifurcation whereby the minimum of the effec- 
tive potential changes stability: it becomes from stable 
center an unstable spiral, with the eigenvalues becoming, 
from purely imaginary, genuinely complex with positive 
real part. Second, if 7 is further increased so as to exceed 
the critical value 



the character of the instability changes via another bi- 
furcation: the fixed point, from unstable spiral becomes 
an unstable node, with the eigenvalues turning, from gen- 
uinely complex, to real positive (via a collision on the real 
axis). However, the global analysis that will be discussed 
in the next section will reveal, in the regime < 7 < 7 cr , 
a third type of bifurcation: a saddle connection breaking, 
which is relative to the behavior of the unstable manifolds 
of the saddle point x . 

To conclude this section, we will now consider the case 
E/2AkkVo = 1, as well as the scenario where no fixed 
points exist. In the first case, there exists solely one fixed 
point, say (5, 0), which occurs from the collision between 
(x_,0) and (x+,0) when a = 0. Then, the Jacobian of 
Eq. (fl~3|) has eigenvalues Ai = and A 2 = (2/3)7/x > 0. 
In this case, the linear stability analysis is insufficient to 
provide any information about the stability of x. Both 
the critical case and that of the absence of fixed points 
will be examined by means of global analysis of the sys- 
tem - see below. 



2. Global analysis 

Let us now study the global phase portrait of the sys- 
tem (llip . We will use two fundamental theorems for 
flows in two-dimensional phase spaces, namely the Hovf 
bifurcation theorem (see, e.g., Theorem 3.4.2 in Ref. [38() 
and the Poincare-Bendixson theorem (see, e.g., Theorem 
1.1.19 in Ref. [Hj]), combined with a nullcline analysis. 
The latter, is an analysis of the direction of the flow 
around the curves consisting of the points on which the 
vector field is parallel (y = 0), or perpendicular (x = 0) 
to the horizontal axis y = |40l |. 

First we consider the case E/2AkkVo < 1, where the 
two fixed points x_ = (x_,0) and x + = (#+,0) exist. 
Since the potential energy U(x) [cf. Eq. (flDjl ] is decreas- 
ing, the maxima (minima) x+ mod(7r//c) [x— mod(7r/fc)] 
are in different and decreasing energy levels, with every 
minimum lying between two maxima. Therefore, in the 
Hamiltonian case, 7 = 0, the conditions for the exis- 
tence of homoclinic orbits (cf. Lemma 14.3 in Ref. (4l| ) 



are satisfied; thus, at the energy levels of the maxima 
x + mod(7r/fc) of U(x), and in the direction of increasing 
energy, correspond to homoclinic orbits on the saddles 
x + mod(7r/fc), surrounding the centers x_ mod(7r/fc). 
Let us denote by W s (5c + ) and W u (5c + ), the stable and 
the unstable manifold of x + , respectively. We may per- 
form a nullcline analysis in the rectangular region Bq = 
{Q < x < f , Y x < y < Y 2 }, where Y 1 < and y = Y 2 > 
are arbitrary horizontal lines, below and above the ho- 
moclinic orbit of x+, respectively. This analysis shows 
that the non-homoclinic branch of W u (5c + ) escapes from 
Bq and becomes unbounded. On the other hand, in the 
Hamiltonian case 7 = 0, the orbits are symmetric with 
respect to the horizontal axis y = 0. The symmetric 
orbit of the unbounded branch of W u (x+) is the non- 
homoclinic branch of W s (x + ) 1 which enters in the region 
Bq, by crossing y = Y2 from below. 

The transition from 7 = to finite < 7 < 7 cr , can 
be analyzed with the help of the Hopf bifurcation theo- 
rem and the Poincare-Bendixson theorem. The former, 
establishes the non-existence of periodic orbits bifurcat- 
ing from the unstable spiral x . This is due to the fact 
that, while the two conditions of the theorem concerning 
the eigenvalues are satisfied, the third one concerning the 
degeneracy condition is not (see, e.g., Ref. 38]). Besides, 
the nullcline analysis in region Bq shows that, apart from 
the stable manifold W s (x+), all orbits escape form the 
strip Y\ < y < I2 , and there is not any other closed orbit 
in the strip. 

Since a periodic orbit does not exist, one may apply the 
Poincare-Bendixson theorem on the invariant manifold 
W s (5c + ). This way, we can establish that each branch 
lying in the half-plane y > joins the unstable spiral 
x+. The situation is more intriguing for the branch of 
W s (x+) lying in the half-plane y < 0, since two possibili- 
ties exist regarding its connection to a fixed point. These 
possibilities depend on the size of < 7 < j CI — 3y/a/fi, 
and on the fact, that the transition from 7 = to fi- 
nite 7 > is a homoclinic bifurcation for the saddle x+, 
where each homoclinic orbit for 7 = breaks. 

When 7 < 7 cr is sufficiently small, a nullcline anal- 
ysis for 7 — > 0, shows that the gradient of the branch 
W s (5t + ) lying in the half-plane y < is close to that 
of its analogue for 7 = 0; furthermore, for 7—^0, the 
spiraling growth is weak. Therefore, since the branch of 
WK s (x+) lying in the half-plane y > is connected to the 
spiral point x_, the other branch of W s (x+) lying in the 
half-plane y < 0, is the one that connects to the saddle 
point x+,1 in the second period of the effective potential 
n/k < x < 2n/k. Hence, for small values of 7 < 7 cr , the 
homoclinic bifurcation creates a saddle connection (see, 
e.g., Refs. [H,|4l|) between x + and x+,1. In this connec- 
tion, the branch of W s (x + ) lying in the half-plane y < 
is - simultaneously - the branch of the unstable manifold 
W^"(x+ ) i) lying also in the half-plane y < 0. 

When 7 < 7 cr is further increased, the branch of 
W s (x+) lying in the half-plane y < 0, connects to 
the unstable spiral point x_ ( i in the second period of 
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the effective potential. Thus, another type of bifurca- 
tion occurs, when 7 < 7 cr exceeds a threshold value 
7thr < 7cr, namely, a saddle connection breaking bifurca- 
tion. The existence of 7thr is guaranteed by the fact that 
as 7 < 7 cr is increased, the growth and oscillation fre- 
quency of the spiraling orbit far from the equilibrium x+ 
dominates, enforcing the considered branch of VF s (x+) to 
form the shape of the spiral. Furthermore, in the regime 
7thr < 7 < 7cr, both branches of W u (5t + ) become un- 
bounded. 

When the bifurcation for 7 > 7 cr occurs, the situation 
is exactly as in the saddle connection breaking setting, 
with the difference that the unstable spirals are replaced 
by the unstable nodes, and the orbits are escaping (con- 
verge) from the strip Y\ < y < Y2 (on the stable mani- 
folds) at an exponential rate. 

We conclude with the case 2 A^kv a ~ ^' wnere solely 
one fixed point x exists, and the case where no- fixed 
points exist. For the former, which is a case of saddle- 
node bifurcation, the nullcline analysis in a small neigh- 
borhood of x reveals that it is unstable. For the latter, 
the nullcline analysis in the region Bo, establishes that 
any orbit must become unbounded. 



3. Phase-plane portraits 

The above analysis on the structure of the phase space 
of the system (fTTj) is illustrated in Fig.[TJ for four different 
values of the damping parameter 7 (other parameters 
values are k = 0.1, E = 0.004, V = 0.05 and fi = 1). 
The computational window is shown for two periods of 
the effective potential U(x). 

In particular, the top left panel presents the family 
of homoclinic orbits in the Hamiltonian case (7 = 0), 
surrounding the stable centers. The non-homoclinic 
branch of W u (x+) diverges and becomes unbounded, 
while W s (x + ) remains bounded. 

In the top right panel, where the orbits have been plot- 
ted for 7 = 0.005, we observe the transition from purely 
imaginary eigenvalues to complex eigenvalues. The sta- 
ble centers of the previous case, now become unstable 
spirals. Furthermore, we see the homoclinic bifurcation 
for the saddle. It is also interesting to observe the nu- 
merically computed orbits in the half-plane y < which 
connect to the saddles, clearly illustrating that the sys- 
tem is in the regime of the saddle connection. 

The bottom left panel corresponds to an increased 
value of 7, namely 7 = 0.05. In this case, the saddle 
connection breaking bifurcation has already taken place. 
Now, the stable manifolds of the saddles consist of orbits 
which are coming from unstable spirals. The unstable 
manifolds, as well as all other orbits, become unbounded 
with increased growth, if compared with the previous 
case. 

Finally, as shown in the bottom right-panel, the solu- 
tions grow exponentially fast in the case 7 = 0.2. The 
character of the instability has changed via the colli- 
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FIG. 1: (Color online) Phase planes for E = 0.004, V = 
0.05, k — 0.1, fi = 1; the computational window is shown 
for two periods of the effective potential U(x). Top left panel 
corresponds to the Hamiltonian case, 7 = 0, while for the 
top right, bottom left and bottom right panels, 7 = 0.005, 
7 = 0.05 and 7 = 0.2, respectively. Bold (red) points depict 
the fixed points, and solid gray (orange) lines show the null- 
clines; the gray (green) lines in the bottom left panel depict 
the homoclinic orbits. 

sion of the two complex eigenvalues which have now be- 
come real, and the unstable spirals have thus become 
unstable nodes. The numerical computed phase plane 
demonstrates the particular node-orbits joining the sad- 
dle points. 

We note that above we have focused on the anti- 
damped case 7 > 0, which is physically relevant in the 
context of BECs. Nevertheless, the dynamical system 
(fTTj) exhibits a rich behavior for 7 < 0, a case correspond- 
ing to a damped Josephson junction (see, e.g., Refs. (42h - 
|44j). As an additional comment along this vein, suitable 
parallels between the two cases can be drawn on the ba- 
sis of the invariance of Eq. ((SJ under the transformation 
7 — > —7 and t — > —t. 

C. Bogoliubov-de Gennes (BdG) analysis 

We complement this dynamical systems picture by 
briefly describing our methods for obtaining the numer- 
ical results for the DGPE (see next section) and com- 
paring them with our analytical predictions. First, we 
will identify the fixed points of the dynamical system 
of Eqs. pT|) - corresponding to stationary dark soli- 
ton states of the DGPE, Eq. fTJ - through a Newton- 
Raphson-type iteration method. Subsequently, we will 
study the linear stability of the stationary solutions of 
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the DGPE by means of the following ansatz: 

u(x, t) = u (x) + e[a(x)e xt + b*{x)e x " t ]. (18) 

Here, uq{x) is the stationary solution of Eq. (TT|) and the 
term proportional to e (which is a formal small param- 
eter) represents a small perturbation. Substituting the 
above expression into Eq. (TTJ), and linearizing with re- 
spect to the amplitudes a and b* , we obtain at 0(e) the 
following linear [Bogoliubov-de Gennes (BdG)] stability 
equations: 



Li -ug 
uf U 



= A 



-2 + 7 






-«-7 



±<9 2 

2 a; 



where Lx = ±c£ — 2\u \ + fi — V and L 2 
2 1 uo | 2 — fx + V. This way, the linear stability problem is 
transformed into a spectral problem for the linearization 
operator: 



\ — i— 7 —2—7 / 

The latter operator is discretized (on the same grid in 
which the solution is computed) and its eigenvalues A 
and eigenvectors {a(x),b(x)} are computed numerically. 

In the following section we will first identify the 
fixed points in the partial differential equation (PDE) 
of Eq. ([T]), and in the (anti-damped Josephson junction) 
equation of motion of the dark soliton center, namely 
the ordinary differential equation (ODE) of Eq. (jSJ). The 
result of the fixed point scheme convergence for the 
PDE equilibrium will be compared to the solution of the 
transcendental equation for the dark soliton center [cf. 
Eq. (|12j) ]. Then, we will compare the eigenvalue prob- 
lem of the linearized ODE [i.e., the explicit expression 
for the eigenvalues of the Jacobian of Eq. (TT3]) ] and that 
of the full DGPE (i.e., the pertinent eigenvalues of the 
operator L). Finally, we will proceed to an even more 
demanding comparison, namely that of the full dynam- 
ics of the infinite degrees of freedom system (the DGPE) 
and that of its one degree of freedom reduction (the ODE 
for the dark soliton center) . Both of these evolutions will 
be followed by means of a time-stepping scheme (namely 
an explicit fourth-order Runge-Kutta method) and the 
results will be illustrated, for direct comparison, in the 
same plot . 



III. NUMERICAL RESULTS 

First, let us offer a few technical remarks, before we 
commence the presentation of our results. We have used 
no flux boundary conditions (as periodic and Dirichlet 
ones seem less appropriate for our setting). Also, in our 
numerical continuations for the exact solutions and their 
linear stability, we have always commenced from a well- 
known, benchmark limit (e.g., the one without potential 
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FIG. 2: (Color online) The top left panel shows the sta- 
tionary solution in (blue) solid line with the dark soliton at 
the local minimum of the external potential shown in (red) 
dotted line; the (green) dashed line shows the background. 
The top right panel shows the corresponding eigenvalue of 
the linear stability problem, with the (blue) star indicating 
the theoretical eigenfrequency (and the (red) circles illustrat- 
ing the DGPE result). The bottom left panel shows the dy- 
namics of the DGPE with the dark soliton initially located at 
xq — —13.89; the (green) dashed line is the evolution of the 
ODE. The bottom right panel shows the DGPE dynamics 
with a dark soliton initially placed at xo = —20.81 (i.e., near 
the turning point for the homoclinic orbit), while the (green) 
dashed line shows the corresponding dynamics of the ODE. 
Parameter values are: 7 = 0, E = 0.004, Vb = 0.05, k = 
0.1, n = \. 



or without the thermal effects) and subsequently per- 
formed parametric continuations in parameters such as 
typically E, Vb and 7. When solutions were found to be 
unstable, the dynamical evolution of the instability was 
monitored through the time dynamics. Our principal re- 
sults can be summarized as follows. 



A. The Hamiltonian case 

1. Solitons located at the center fixed point 

In the Hamiltonian case, i.e., for 7 = 0, the phase 
plane of Eq. ([TT]) is shown in the top panel of Fig. Q] 
(other parameter values are: E = 0.004, Vb = 0.05, k — 
0.1, fi = 1). Direct inspection of the phase plane sug- 
gests that the center of the dark soliton may periodi- 
cally oscillate around the (center) fixed point (x_,0) = 
(—13.6363,0) with a linearization frequency (of small 
amplitude oscillations) of u) e g = y/a — 0.0302 (re- 
call that a = kyjAAj,k 2 V£ — E 2 ). On the other hand, 
there is a homoclinic orbit passing through the saddle 
(x+,0) = (—2.0709,0) which intersects the x-axis at (a 
finite turning point of) about —20.5. This phenomenol- 
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FIG. 3: (Color online) The top left panel shows the stationary 
solution with the dark soliton at the maximum of the effec- 
tive potential. The top right panel shows the corresponding 
eigenvalues of the linear stability analysis; Again, in addition 
to the (red) circles stemming from the eigenvalues of the op- 
erator L, there exists the theoretically predicted growth rate 
[indicated by the (blue) star] from the corresponding ODE. 
The two bottom panels show the dynamics of the DGPE, with 
the dark soliton initially placed at the local maximum of the 
effective potential, with two different perturbations: the dif- 
ference refers to the sign of the unstable eigenvector of the 
linearization problem, which has been added to the unstable 
dark soliton in order to initiate the instability evolution. The 
amplitude used in the two cases shown is ±0.1. Parameter 
values are: 7 = 0, E = 0.004, V = 0.05, k = 0.1, fx = 1. 



ogy at the level of the DGPE (and the comparison with 
the ODE) is shown in Fig. [3J In the top left panel of 
the figure, the solid (blue) line shows an example of the 
exact stationary solution at the effective potential mini- 
mum as obtained from our fixed-point iteration method. 
Furthermore, the dotted (red) line denotes the external 
potential, while the dashed (green) line is the the back- 
ground Ub(x), i.e., the ground state of the system. The 
latter showcases the validity of our ansatz approxima- 
tion assuming that the dark soliton can be embedded in 
the system through a product ansatz, where it is directly 
multiplied to the (back)ground state of the system. The 
top right panel of the figure presents a comparison of 
the eigenvalues of the matrix L which are found to be 
purely imaginary, corroborating the neutral stability of 
this center point. Within this plot, the (blue) star rep- 
resents the theoretical frequency ±0.0302, which nearly 
coincides with one of the eigenvalues ±0.0292i of L. As 
has been analyzed in detail in previous works (see, e.g., 
Refs. [2|, [ll|, |45|) for the Hamiltonian case, and even for 
the thermal BEC setting [2(J this eigenvalue rep- 
resents the so-called anomalous (or negative energy, or 
negative Krein signature) @ mode, which is the one that 
pertains to the (eigenfrequency of the) motion of the dark 
soliton around its equilibrium position. As seen in the 



FIG. 4: (Color online) The top and bottom sets of panels are 
for the dark soliton at the local minimum and maximum of the 
effective potential, respectively. The left panels showcase the 
change of frequency (or unstable eigenvalue) vs. Vo for E = 
0.004, while the right ones illustrate the change of frequency 
(or unstable eigenvalue) vs. E for Vo = 0.05. The (red) dashed 
line is the result obtained analytically from the ODE, while 
the solid (blue) line represents the corresponding mode within 
the spectrum of linearization around the dark soliton of the 
DGPE. Parameter values are: 7 = 0, k = 0.1, y, = 1. 

bottom left panel of Fig. [21 the (white) dotted line os- 
cillating with the theoretical frequency 0.0292 provides 
a quite accurate description of the full dynamics of the 
DGPE around its center equilibrium position (the soliton 
is initially placed at xq — —13.89), confirming our expec- 
tation that it is indeed this eigenmode that describes the 
dark soliton motion. Nevertheless, we have subjected the 
ODE approximation to a far more stringent test by ini- 
tializing the soliton at xq = —20.81, namely very near 
the turning point for the homoclinic orbit (i.e., the most 
unstable orbit of the system). However, even in that case 
as attested by the bottom right panel of the figure, the 
agreement between the ODE and PDE results is quite 
striking. 

2. Solitons located at the saddle fixed point 

We now turn to an examination of the unstable saddle 
fixed point existing at the local maximum of the effective 
potential, as illustrated in Fig. [31 The relevant solution 
is again shown in the top left panel, while, this time, the 
linear stability problem of the DGPE possesses a real pair 
of eigenvalues illustrating its saddle character. These are 
again in reasonable agreement with the expectation of 
the perturbation theory growth rates, as shown in the 
top right panel of the figure. In order to show differ- 
ent possibilities of evolution associated with this saddle 
point, we perturb its unstable eigenvector with a positive 
amplitude and with a negative amplitude, respectively, 
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in the two bottom panels of Fig. |3l In the case of the 
bottom left, the dark soliton acquires a velocity towards 
positive values of x and thus escapes to infinity sliding 
downhill (i.e., washboarding) through the subsequent po- 
tential wells. In the case of the right panel, the opposite 
sign of the perturbation forces the soliton to move in the 
opposite direction, where it encounters a turning point. 
However, after hitting the turning point and changing 
direction of motion, it has enough energy to overcome 
the local energy maximum (at the saddle point) and also 
washboard towards +00. In both cases, we find that 
the ODE dynamics captures fairly accurate the observa- 
tions of the DGPE. Nevertheless, in the right one the 
dynamics is very close to the homoclinic orbit and hence 
is most unstable, a feature which is observed with the 
longer meandering of the ODE trajectory in the vicinity 
of the saddle point before it eventually escapes towards 
+00. It should be noted that this feature should perhaps 
naturally be expected: in the infinite degrees of freedom 
system (the PDE), the soliton has the ability to shed 
(weak, but important in this case) radiative wavepack- 
ets, whereby it becomes shallower and faster and, hence, 
it is easier for it to overcome the relevant barrier. On the 
other hand, such decay channels are forbidden by virtue 
of the ansatz in the one degree of freedom reduction. 

Having captured an individual example of both types 
of fixed points, we illustrate in Fig. @] the results of a 
corresponding continuation over the two parameters, the 
strength E of the linear (electric or gravitational) poten- 
tial and that of the periodic (OL) potential, Vq, in both 
cases for 7 = i.e., in the Hamiltonian case. As pre- 
dicted by the perturbation theory, when the dark soli- 
ton is around the local min (max) of the effective poten- 
tial, the frequency (growth rate) of the movement of dark 
soliton increases as Vq increases, while it decreases as E 
increases. Both of these dependencies are essentially en- 
capsulated in the functional form of a (since the relevant 
eigenvalue is predicted to be ±ia in the top panel, while 
should be ±a in the bottom one. 



3. On the accuracy of perturbation theory 

It is well-known that in the case of the periodic (at 
infinity) linearization operator that emerges in the pres- 
ence of solely an optical lattice potential, the spectrum 
acquires bands separated by gaps. On the other hand, 
the linear potential due to its unbounded nature at ±00 
renders the linearization spectrum of eigenstates again a 
purely point spectrum (as in the more customary case 
of a parabolic trap). The above discussion justifies the 
pure point spectrum nature of the excitation frequencies 
observed in Figs. [5] and [3] However, as is well-known vari- 
ations of the parameters of the system (e.g., as was seen 
above, the increase of Vq in the stable center case of the 
effective potential minimum) can have the consequence 
that they increase the soliton's linearization frequency, 
which corresponds to its oscillation frequency inside the 




FIG. 5: (Color online) The top left panel shows the station- 
ary dark soliton located at the local minimum of the effec- 
tive potential. The top right panel shows its corresponding 
eigenvalues of the linear stability problem, which indicate 
its oscillatory instability. The bottom panel shows the full 
dynamics of Eq. (fT| in its underlying contour plot with the 
dark soliton initially at xo = —3.41 (around the center fixed 
point). The dynamics of the ODE for the same setting is 
indicated by the dashed (green) line. Parameter values are: 
7 = 0, E = 0.03, Vo = 0.1, k = 0.5, (jl = 1. 



shallow trap of its local well in Fig. [2] As this parame- 
ter increases, there will emerge critical points where the 
anomalous mode (of negative energy) associated with the 
dark soliton will collide with the positive energy modes of 
the background. These collisional events in the spectrum 
will produce instabilities 0, EH of the oscillatory type as- 
sociated with complex eigenvalue quartets (when 7 = 0). 
An example of this type is observed in Fig. O This is 
a feature of the full DGPE problem that our reduction 
to a one degree of freedom ODE cannot capture. Hence, 
we highlight this potential deficiency of the theory, which 
perceives the soliton as a stable entity, while the dynam- 
ics of the bottom panel of the figure clearly illustrate the 
coherent structure to be oscillatorily unstable and hence 
eventually departing from its local well. 

Another parameter variation that renders the pertur- 
bation theory less accurate concerns the increase of k. 
This feature is illustrated in the top right panel of Fig. [51 
for k = 0.5 and the dark soliton being located at the 
local maximum of the external potential. The difference 
between the unstable eigenvalue stemming from the lin- 
earization of the PDE and that from the ODE is obvi- 
ously nontrivial, and is also clearly mirrored on the bot- 
tom panel representing the dynamics of the ODE and the 
DGPE. The over-estimate (by the ODE) of the growth 
rate of the unstable eigenmode leads the dark soliton 
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FIG. 6: (Color online) The top left panel represents the 
stationary dark soliton located at local maximum of effec- 
tive external potential. The top right panel shows the cor- 
responding eigenvalues of the linear stability problem via 
(red) circles, while the (blue) star represents the less suc- 
cessful in this case analytical prediction. Finally, the bot- 
tom panel shows the dynamics of Eq. |T} with a pertur- 
bation of the unstable saddle point through its eigendirec- 
tion of amplitude of 0.1; The dynamics of the ODE is in- 
dicated by the dashed (green) line. Parameter values are: 
7 = 0, E = 0.03, Vo = 0.1, k = 0.5, fj, = 1. 



"particle" to diverge from the unstable equilibrium point 
much faster in the particle approach than in the full PDE 
one. This can presumably be attributed to the emerging 
competition between the length scale of the lattice and 
that of the soliton proper (see discussion in Ref. HU), 
which contributes towards flaws of the dark soliton per- 
turbation theory for moderate values of k. 

From now on, we will restrict our considerations to the 
regime of small k (in particular, fixing k = 0.1), so as to 
ensure the validity of the perturbation theory and avoid 
the competition of the soliton and lattice length scales. 



B. Finite-temperature-induced dynamics 

We will now turn to a detailed examination of the role 
the temperature-dependent parameter 7. A direct obser- 
vation concerning our dynamical system in the latter case 
is that the Hamiltonian %(x,y) = \y 2 + U (x), instead of 
being conserved as in the case, 7 = 0, rather satisfies 



m 

dt 



dy 
dt 



y - 



dx 



y 



> 0. 



(19) 



Thus, the phase portrait instead of consisting of the level 
sets of a constant energy U (as in the case of 7 = 0), will 
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FIG. 7: (Color online) The case of the dark soliton at local 
minimum of the effective potential. Three case examples of 
the spectral plane of the imaginary versus real part of the 
relevant eigenvalue A = Re(A) + ilm(A) is shown. The top 
left panel is for 7 = 0.005 < 7 cr , the top right panel is for 
7 = 0.0905 = 7cr, and the bottom right panel is for 7 = 0.1 > 
7c r ; (blue) stars depict the analytically derived eigenvalues 
of Eq. p5p. Other parameter values are: E — 0.004, Vo = 
0.05, k = 0.1, fi = l. 



rather consist of curves that tend to higher levels of U 
during the time evolution. 



1. Sohtons located at the local minimum 

First, we consider the case of a stationary dark soliton 
located at the minimum of the effective potential. As 
indicated above (sec. II. B), with an increase of 7, the 
stationary solutions of the ODE and the DGPE remain 
unchanged, yet the linear stability properties of such a 
solution change upon crossing the critical value 7 cr [cf. 
Eq. (|17p]: the local minimum of the effective potential 
becomes, from oscillatorily unstable (for 7 < 7 cr ), expo- 
nentially unstable (for 7 > 7 cr ); these changes of stability 
are mirrored in the phase plane portraits of Fig. [T] 

First, employing the BdG analysis, we will investigate 
if the above ODE-based stability picture complies with 
the stability of the (formerly) anomalous mode of the 
soliton placed at the local minimum of the effective po- 
tential. In our simulations, we use the parameter values 
E = 0.004, V = 0.05, k = 0.1, fi = 1, resulting in 
the value 7 cr = 0.0905. For these values, in Fig. [7] we 
illustrate the spectral planes obtained by the BdG anal- 
ysis, for three different values of 7 (below, at, and above 
7 cr ). It is clear that, due to the dissipative type of dy- 
namics of the DGPE, all eigenvalues of the operator L 
tend to the left half of the complex plane (thus having 
a negative real part and corresponding to decaying ex- 
citations) except for the anomalous mode. The latter, 
as justified in Ref. [46], has the opposite behavior and 
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tends to right half of the complex plane illustrating the 
unstable dynamics of the excitation. In the figure, it is 
observed that the analytically predicted relevant eigen- 
values of Eq. (|15p [(blue) stars] are in excellent agreement 
with the anomalous mode eigenvalues [(red) circles] ob- 
tained numerically via the BdG analysis. 

As an additional interesting observation, we note the 
observable bending of the stable modes in the left half of 
the complex spectral plane in Fig. [7] and the formation of 
two curves of eigenvalues for larger values of 7 presum- 
ably associated with the two different asymptotic states 
within our linear plus periodic potential. 

The above mentioned feature of the anomalous mode 
(i.e., its motion toward the right half plane) is deeply con- 
nected to the fundamental difference between dynamic 
and thermodynamic instability of the dark soliton: the 
soliton is an excited state of the condensate and, hence, 
while dynamically stable for the Hamiltonian variant of 
the problem, it is thermodynamically unstable. However, 
in order to realize this instability and enable the decay of 
the excited system to its ground state, the role of dissi- 
pation (induced by the presence of a finite, temperature- 
dependent parameter 7) is essential. This is mirrored in 
the direct instability of the anomalous mode, as shown in 
Fig. [5J in excellent agreement with our ODE predictions 
for both the real and the imaginary part of the relevant 
complex eigenvalue. Notice also that even the critical 
point for the collision of the complex pair of eigenvalues 
and its subsequent transformation into a pair of positive 
real eigenvalues (when the fixed point transforms from an 
unstable spiral into an unstable node) is very accurately 
captured. 
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FIG. 8: (Color online) The figure shows the real (top panel) 
and imaginary (bottom panel) part of the unstable eigenvalue 
pertaining to the dark soliton located at the effective potential 
minimum vs. 7. The (red) dashed line is the result of the 
perturbation theory analysis, while the (blue) solid line is 
the corresponding mode of the DGPE. Parameter values are: 
E = 0.004, Vo = 0.05, k = 0.1, fx = 1. 



2. Solitons located at the local maximum 

We now turn to the examination of the effect of tem- 
perature in the case of a soliton located at the local max- 
imum of the effective potential. As predicted also by the 
ODE, the linear stability problem for the PDE in this 
case has one unstable eigenvalue (i.e., the relevant fixed 
point is once again a saddle). As can be observed in 
Fig. |9l the agreement is once again good for sufficiently 
small values of 7, but becomes less adequate as 7 in- 
creases. Nevertheless, as can be observed, e.g., in Fig. 11 
of Ref. 22], physically relevant values of (the dimension- 
less) parameter 7 are of 0(1O -3 ) for realistic BECs, hence 
our approximation is quite reasonable in that range. 

Finally, as seen in Fig.[T0l the dynamics of the unstable 
dark soliton in the presence of the effect of temperature 
(and for values of 7 in the physically relevant range as 
discussed above) is in very good agreement between the 
DGPE and the corresponding ODE. This is illustrated 
both in the oscillatorily unstable, for this range of 7's, 
case of the minimum of the effective potential and in the 
exponentially unstable case of its maximum. 
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FIG. 9: (Color online) The figure shows the stability charac- 
teristics of the dark soliton for the case of the local max of the 
effective potential. The left panel presents a typical example 
of the eigenvalues of the spectral plane for 7 = 0.005. The 
right panel shows the real unstable eigenvalue of the relevant 
saddle point vs. 7, with the (red) dashed line illustrating the 
theoretical prediction, while the solid (blue) line is the corre- 
sponding linear stability eigenmode of the DGPE. Parameter 
values are: E = 0.004, Vo = 0.05, k = 0.1, fj, = 1. 



IV. CONCLUSIONS 

We have studied matter-wave dark solitons in atomic 
Bose-Einstein condensates with a linear and periodic (op- 
tical lattice) potential in the absence as well as in the 
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FIG. 10: (Color online) The figure presents contour plots 
of the dynamics of the DGPE for parameter values 7 = 
0.005, E = 0.004, Vo = 0.05, k = 0.1, fj, = 1. The left 
panel shows the case of the dark soliton initialized at -14.77 
in the vicinity of the effective potential energy minimum, with 
the (green) dashed line showcasing the results of the pertur- 
bation theory. The right panel is for the dark soliton initially 
located in the vicinity of the maximum of the effective poten- 
tial (at —0.414). Here once again, the (green) dashed line is 
from the evolution of our perturbative ODE. 

presence of thermal effects. Our fundamental model 
employed in the analysis was the dissipative Gross- 
Pitaevskii equation. Using dark soliton perturbation the- 
ory, we found that the motion of the dark soliton cen- 
ter is governed by an ODE, namely a Newtonian equa- 
tion of motion, with an anti-damping term; this ODE is 
strongly reminiscent of the one describing the dynamics 
of the phase in a superconducting Josephson junction. 
We found that for sufficiently small wavenumbers of the 
periodic potential, such that the length scale of the soli- 
tary wave is much shorter than that of the lattice, the 
ODE gives a very accurate prediction for the evolution 
of the center of the dark soliton. On the other hand, for 
larger wave numbers (where a competition of the scale 
of the soliton and that of the periodic potential would 
emerge) , or for sufficiently strong potentials (beyond the 
realm of perturbative effects), naturally, the approach 
was found to be less accurate. 

In the case of Hamiltonian dynamics (zero- 
temperature limit), the resulting effective potential 



was found to possess effective minima that were typi- 
cally stable - although potentially oscillatorily unstable 
- and effective maxima, corresponding to saddle points. 
The saddle points preserved their nature under thermal 
perturbations, while the centers became unstable spirals 
and, finally - for sufficiently large strength of relevant 
perturbation - unstable nodes. These features (with the 
exception of the oscillatory instability of the Hamiltonian 
case that requires a higher-dimensional description to 
be captured) were all adequately represented by our 
equation of motion. In addition to the statics and 
near-equilibrium linearization, the ODE was found to 
yield a very good approximation of the full PDE in the 
corresponding dynamics. 

There are many directions that one can envision as 
potential themes for future study. On the one hand, a 
detailed examination of the ground state of the system 
and especially of its linearization spectrum would be a 
first step towards considering the "interaction" of this 
spectrum with the anomalous mode of the solitary wave. 
To address such an interaction, an expanded higher num- 
ber of mode ansatz is necessary in order to capture the 
relevant resonance effects. Another direction that can be 
envisioned concerns a consideration of dynamical systems 
involving two or more solitons, in the spirit of Ref. [ll| in 
the absence, and [HJ in the presence of thermal effects. 
For more such solitary waves, more complex phenomena 
including the possibility of chaotic orbits may become 
manifest. Another generalization may concern the study 
of two-dimensional such settings, where the study of vor- 
tices and their evolution would be of interest. These top- 
ics will be presented in future works. 
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